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ABSTRACT 

I 

The  governing  equations  for  constrained  multibody  systems 

.ire  formulated  in  i  wunne r  suitable  their  automated ,  nu-\.  ,-i.  .il 
development  and  solution.  Specif icaily ,  the  "closed  loon"  proniem 
of  multibocy  chain  systems  is  addressed. 

The  governing  equations  are  developed  by  modifying  dynamical 
equations  obtained  from  Lagrange's  form  of  d'Alembert's  principle. 

T’ais  modification,  which  is  based  upon  a  solution  of  the  constraint 
equations  obtained  through  a  "zero  eigenvalues  theorem,"  is,  in 
effect,  a  contraction  of  the  dynamical  equations. 

It  is  observed  that,  for  a  system  with  n  generalized  coordinates 
and  m  constraint  equations,  the  coefficients  in  the  constraint 
equations  may  be  viewed  as  "constraint  vectors”  in  n-dimensional 
space,  Then,  in  this  setting  the  system  itself  is  free  to  move  in 


the  n-m  directions  which  are  "orthogonal"  to  the  constraint  vectors. 
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INTRODUCTION 


This  report  presents  a  formulation  of  the  governing  equations  of 
constrained  multibody  systems.  The  objective  is  the  establishment  of 
procedures  for  the  automated  generation  of  the  equations. 

Recently  there  has  been  an  increasing  interest  in  the  efficient 
development  of  governing  dynamical  equations  of  multibody  systems.  This 
interest  is  stimulated  by  the  fact  that  manv  physical  s vs  tens  can  be 
nio .  ...  j  :v.r  sysCcfi’s  i  >.  r.r.oct^c  n . "  j  •  .'vi  itv  .  '  *■  -r:.’  'C  .t.'.'w 

physical  svstems  of  interest  are  robots,  manipulators,  human  bodv  models 
and  biodynamic  systems,  rad  flexible  cables  or  '-bains. 

There  have  been  a  numoer  of  formulations  of  the  equations  of  motion 
of  multibody  systems  [1-19]*.  The  majority  of  these  have  been  restricted 
to  "open  chain"  or  "open  tree"  systems:  that  is,  systems  of  rigid  bodies 
such  that  adjacent  bodies  have  at  least  one  common  point  and  such  that 
no  closed  loops  are  formed.  Figure  1.  illustrates  such  a  system.  The 
formulation  of  the  governing  equations  of  motion  of  such  systems  has 
advanced  to  the  point  where  the  coefficients  of  the  governing  differential 
equations  can  be  formed  automatically  (numerically)  by  simply  knowing  the 
connection  configuration  [10-12] . 

However,  during  recent  years,  there  has  also  been  interest  in  the 
dynamics  of  systems  possessing  closed  loops,  where  some  of  the  branches 
of  the  tree  or  chain  are  connected.  Figure  2.  illustrates  such  a  system. 
These  svstems  are  useful  in  modelling  such  physical  systems  as:  closed 

*  Numbers  in  brackets  refer  to  References  at  the  end  of  the  report. 
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mechanisms,  "docking"  manipulators  of  spacecraft,  snip  cranes,  restrained 
human  body  models,  and  cables  anchored  at  both  ends. 

As  noted  above,  this  report  presents  a  procedure  for  the  automatic 
formulation  of  the  governing  equations  of  such  closed- loop  multibcdy 
systems.  The  procedure  is  based  upon  Lagrange's  form  of  d'Alembert's 
principle  as  exposited  by  Kane  et.  al.  [14,20-22]  and  as  used  in  [9-12] 
to  develop  the  dynamical  equations  of  motion.  It  is  also  based  upon  a 
"zero  eigenvalues  theorem"  as  exposited  by  Walton  and  Steeves  [23]  to 
provide  an  automatic  inclusion  in  the  analysis  of  the  constraint  equations 
The  balance  of  the  report  itself  is  divided  into  parts  with  the 

following  part  providing  some  preliminary  information  useful  in  the  sequel 
This  includes  a  review  of  dynamical  formulations  of  multibody  systems  and 
a  statement  of  the  "zero  eigenvalues  theorem."  This  is  followed  in  the 
next  part  by  the  governing  equation  formulation  for  constrained  or  closed- 
loop  multibody  systems.  The  subsequent  part  presents  a  simple  example. 

The  final  two  parts  discuss  generalizations  and  other  features  of  the 
formulation. 

PRELIMINARY  CONSIDERATIONS 

Coordinates  and  Kinematics 

Consider  again  the  multibody  system  of  Figure  1.  This  system  will 
have,  in  general,  3N+3  degrees  of  freedom  where  N  is  the  number  of  bodies 
of  the  system.  These  degrees  of  freedom  might  be  delineated  as  follows: 
Arbitrarily  select  a  body  of  the  system  as  a  reference  body.  Call  this 
body  B^.  Next,  label  or  number  the  remaining  bodies  of  the  system  in 
ascending  progression  away  from  B^  through  the  branches  of  the  tree 


structure,  moving  clockwise  from  branch  to  branch.  Then  the  orientation 
of  B,  relative  to  a  fixed  (inertial)  reference  frame  R  together  with  the 
orientation  of  the  remaining  bodies  of  the  system  relative  to  their 
adjacent  lower-numbered  bodies  defines  3N  degrees  of  freedom.  Finally, 
the  location  of  an  arbitrary  reference  point  in  relative  to  R  defines 
an  additional  3  degrees  of  freedom. 

The  position  and  configuration  of  the  system  can  thus  be  described  by 
3N+3  generalized  coordinates  .  Let  y^  (1=1, . . . ,3N+3)  represent  their 
time  derivatives*.  Next,  let  n^  (i=l,2,3'  represent  a  mutually  perpendicular 
.nit  vector  set  fixed  in  R.  _et  G,  represent  the  mass  center  of  oouv  B. 
(k=l,...,N>.  Then,  it  has  been  shown  r 20 , 21]  that  the  velocitv  of  G,  in  R 

rC 

and  the  angular  velocity  of  B  in  R  may  be  expressed  in  the  form: 


V,  =  V  v  n 
~k  kim  1-m 


and 


=  o 


y  n 

kim  I'm 


(1) 


(Regarding  notation,  a  repeated  index,  such  as  l  or  m  in  Equation  (1) 

represents  a  sum  over  the  range  of  that  index,  unless  otherwise  stated.) 

The  coefficients  v,  ,  and  ,  in  Equation  (1),  and  their  derivatives, 
kim  klm 

play  a  central  role  in  the  analysis  of  the  sequel.  They  are  components  of 
the  so-called  "partial  velocity"  and  "partial  angular  velocity"  vectors: 
3vk/9yi  and  3^^/3yr.  These  vectors  are  useful  in  forming  the  generalized 


*  The  reason  for  using  the  symbol  y?  instead  of  x^  is  that  there  exist 
instances  when  a  convenient  choice  of  generalized  coordinate  derivatives 
result  in  functions  vj  which  cannot  be  integrated  to  obtain  the  coordinate 

x  .  In  such  cases,  the  x  do  not,  in  general,  exist  (and  are  sometimes 

t  z 

called  "quasi-coordinates").  This  occurs,  for  example,  when  the  y,  are 

K, 

selected  as  angular  velocity  components.  See  [24] . 
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accive  and  inertia  forces  of  the  system.  The  coefficients  v,  ,  and  ut.  . 

k-.m  klm 

and  their  derivatives  may  be  formed  by  simple  multiplication  algorithms 
as  developed  in  [9-12].  Hence,  by  differentiating  in  Equation  (1),  the 
acceleration  of  G  in  R  and  the  angular  acceleration  of  B,  in  R  may  be 

K.  K 

expressed  as: 


=  (v  y 


+  v.  „  y„)n  and  a,  =  (w,  ,  y„  +  ui,  .  y„)n  (2) 
klnr  J,  ~m  -k  kim  l  k£m  l  ~m 


Equations  of  Motion 

Constcer  the  system  in  Figure  ±.  ■  to  oe  subjecteu  to  an  eucernaa^y 
applied  force  field  which  may  be  represented  on  a  typical  body  B,^  by  a  single 
force  F^  passing  through  together  with  a  couple  with  torque  M^.  Similarly, 
let  the  inertia  force  system  on  be  represented  by  a  single  force  F* 
passing  through  G^  together  with  a  couple  with  torque  M*.  Then  F*  and  M* 
may  be  expressed  as  [21] : 


-m,  a 
k.  k 


(no  sum) 


(3) 


and 


$ 


■h 


?k  ~  ik x 


.o,  )  (no  sum) 
k 


(4) 


where  is  the  mass  of  and  1^  is  the  inertia  dyadic  of  relative  to 
G^.  Through  use  of  orthogonal  transformation  matrices  [10]  1^  may  be 
expressed  in  the  form: 


I.  =  I,  n  n 
-k  kmn~m~n 


(5) 


Lagrange's  form  of  d'Alembert's  principle  then  leads  to  governing 
dynamical  equations  of  motion  of  the  form  [21]  : 
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F .  +  F* 


=0  l  =  1,. . . ,3N+3 


(6) 


wnere  r .  is  ca^ie- 


"generalized  active  "’orce"  and  may  \-e  exp ressed  as: 


F  =  v  F  +  u  m 
t  k£m  km  k?.m  tan 


(7) 


where  there  is  a  sum  from  1  to  N  on  k  and  from  1  to  3  on  m,  and  where 
F^m  and  are  the  components  of  F^  and  Similarly,  F*,  in  Equation 

(6),  is  called  the  "generalized  inertia  force"  and  may  be  expressed  as: 


F*  =  v.  ,_F*  +  .  M*  ($) 

•<. .  .T  v.r. 

where  there  is  a  sum  from  1  to  N  on  k  and  from  1  to  3  on  - ,  and  where 

F?,;  and  M*  are  the  n  components  of  F*  and  M*. 

xm  •  -k  -K 

Substituting  from  Equations  -l'  to  ■'?)  into  (7)  and  »'  •  '•  and 

finally  into  (6)  leads  to  the  equations  of  motion  which  may  be  written 
in  the  form [10] : 


■  q-  q 


l  =  1, . . .  ,3i>3 


<?) 


where  there  is  a  sum  from  1  to  3N+3  on  q  and  where  a,  and  f,  are  given  bv : 

M  q 


a,.q  ^V^k.'m'kqm  ^kmn  k^m  kqm 


and 


f  =  F  -  (m.  v,  v.  v  +  I  o,  ,  o, 
k  km  kura  a  kmn  k.n  k 


v 

un'  u 


6nmh^kmr  J'kun  "ksr^k.  n’’  u-'  s 


un 


where  there  is  a  sum  from  1  to  N  on  k,  from  1  to  3N+3  on  u  and  s  and  from 


1  to  3  on  the  other  repeated  indices  and  where  e  , 

r  nmh 

tion  symbol  [25]. 


is  the  standard  permuta- 


Constraint  Equations 

Equations  (91  represent  the  general  governing  dynamical  equations 
for  open  chain  or  open  tree  systems.  However,  if  the  system  has  one  or 
more  closed  loops,  as  illustrated  in  Figure  2.,  there  are  additional 
equations  which  need  to  be  satisfied  to  insure  that  the  closed  loops  are 
maintained  throughout  the  motion  of  the  system.  These  equations  are 
holonomic  constraint  equations  [21]  and  they  may  be  written  in  the  form: 

g.(x,)  =  0  i  =  l,...,m  :  m  <  3M+3  (12) 

(These  equations  may  be  obtained  by  simply  adding  to  zero  the  relative 
position  vectors  of  the  connecting  joints  around  the  respective  loops.) 

It  should  be  noted  that  constraint  equations  of  the  form  of  Equation  (12) 

can  arise  in  ways  different  than  that  of  the  closed  loops  mentioned  above. 
This  can  occur,  for  example,  with  restrictions  on  the  motion  at  a  joint  or 
with  the  anchoring  of  one  or  several  of  the  bodies  to  a  fixed  frame  R. 
Finally,  by  differentiating,  Equation  (12)  becomes  a  linear  relation  in  the 
v,  and  may  be  expressed  in  the  form: 

by,  =  0  i  =  l,...,m;  i  =  l,...,3N+3  (13) 

where  the  b  are,  in  general,  functions  of  x,  and  t.  Equations  (9)  and 
(13)  thus  constitute  the  governing  equations  for  a  "closed-loop"  system. 
These  are  co  be  cast  into  a  solvable  form  in  the  sequel. 

Zero  Eigenvalues  Theorem 

For  a  constrained  N  body  chain  system,  the  n  dynamical  equations  (9) 
together  with  the  m  constraint  equations  (13)  constitute  n+m  equations  for 
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the  n  unknown  y.,  where  n  is  3N+i.  Hence,  the  system  is  over  determined. 

One  approach  to  overcoming  this  difficulty  is  to  solve  Equations  (\13) 

for  m,  say  the  last  m,  of  the  y,  in  terms  of  the  first  n-m  v,  as  "independent' 

C 

generalized  coordinate  derivatives,  the  partial  velocities  and  the  partial 
angular  velocities  can  be  expressed  exclusively  in  terms  of  these  v.  [21]. 
Finally,  by  following  the  procedure  suggested  by  Equations  (6),  (7),  and 
(8),  n-m  governing  dynamical  equations  are  obtained  for  the  y„  (1=1 , . . ,n-m) . 

Although  this  approach  is  suitable  for  relatively  small  systems, 
are  js  m  l o l  ' .  .  .  line,  it  '  ■;*  sy  s r e".s  —  ~ •  i ..  ^  Lir  Ly  . n 

attempting  to  automate  it.  Among  these  di: .Lenities  is  the  problem  of 
obtaining  a  consistent  solution  ■  Equations  (13)  for  it  of  the  y.  in  terms 
of  the  remain in :  n-m  .  Another  iifficujtv  is  the  problem  of  automatically 
eliminating  these  m  y,  from  the  partial  velocities  and  the  partial  angular 
velocities.  However,  in  1966,  while  working  on  a  constraint  problem  of  a 
different  context,  Walton  and  Steeves  [23]  developed  an  automated  procedure 
for  solving  equations  such  as  Equations  (13),  for  m  of  the  y,  in  terms 
of  the  remaining  independent  n-m  y,.  An  extension  of  their  procedure  can 
be  developed  to  automatically  eliminate  m  of  the  v,  from  the  martial  velocitv 
and  partial  angular  velocity  vectors.  Their  procedure  and  its  extension  are 
based  on  a  "zero  eigenvalues  theorem"  as  outlined  in  the  following  oaragraphs 

Consider  Equations  (13)  to  be  written  in  the  matrix  form  as: 

Bv  =  0  (14) 

where  3  is  an  mxn  rectangular  matrix  with  elements  b .  .  and  v  is  an  n 

i. 
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element  column  matrix  with  elements  v..  From  B  form  the  nxn  symmetric 
matrix  S  defined  ns: 


S  =  B 


(15' 


where  S1  is  the  transpose  of  B.  Since  S  is  symmetric,  there  exists  an 
orthogonal  matrix  T  such  that: 

TTST  =  A  (16) 

where  ..  a  a  nxn  Jiaconai  matrix  with  real  elements  or  "eigenvalues1' 

‘  .  (i  =  l,...,n)  [26],  .'h.ese  eigenvalues  are  readily  seen  to  be  non-negative 

T  t  7 

as  rollows:  Let  v  be  a  typical  column  of  1  and  let  w  be  Bv.  Then  w  w  *  v'B  8v  = 

T  7  ,  .  T 

"  <v.  But  w  _  0,  and  by  Equation  (16),  v  Bv  is  seen  to  be  an  element  of 

A,  say  \ ^ .  Hence,  \  _>  0.  It  is  also  readily  seen  that  there  exist  zero 

eigenvalues:  Since  B  is  an  m x n  matrix,  its  rank  is  less  than  or  equal 
to  m  (26] .  Then,  bv  Equation  (15),  the  rank  of  S  is  also  less  than  or  equal 
to  m.  But,  since  m  <"  n  the  rank  of  the  nxn  matrix  S  is  less  than  n. 

Let  the  columns  of  T  in  Equation  (16)  be  arranged  so  that  the  eigen¬ 
values  of  S,  or  the  diagonal  elements  of  are  ordered.  That  is.  arranse 

T  such  that  '>  >  '•  >  ...  '■  '•  .  (From  the  preceeding  argument,  the  last  r> 

1  ~~  ~  ~  —  n 

of  these  will  be  zero,  where  p  ^  n-m.)  Next,  let  the  mxn  matrix  D  be 
defined  as: 

D  =  BT  (17) 

Then,  from  Equations  (15)  and  (16)  it  is  seen  that: 

DTD  =  A  (13) 


i 

5 

3 

i 


1 
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Hence,  since  Che  last  p  rows  (and  columns)  of  A  are  zero,  C  may  be  written 
i.n  the  partitioned  form: 

D  =  [D| 0]  (19) 

where  D  is  an  (n-p)x(n-p)  matrix  with  mutually  orthogonal  columns,  and 

T 

where  n-m  <_  p  <  n.  By  noting  that  TT  is  the  nxn  identity  matrix,  the 
constraint  equation  (14)  may  be  written  as: 

Dy  =  DTTTy  2  ;yrz  =0  (20) 

where  z  is  the  n  element  column  matrix  defined  as  Tiv.  In  view  of  Equation 
(19),  the  final  equality  in  Equation  (20)  is  satisfied  if  the  first  n-p 
elements  of  z  are  zero,  irrespective  of  the  values  o:  the  last  p  elements 
of  z. 

Since  T  is  orthogonal,  the  definition  in  Equation  (20)  may  be 
"inverted"  leading  to  the  expression: 

y  =  Tz  (21) 

However,  since  the  first  n-p  elements  of  z  are  zero,  y  may  be  rewritten  as: 

y  =  Tz  (22) 

where  T  is  the  nxp  matrix  whose  columns  are  the  last  p  columns  of  T. 

(In  view  of  the  ordering  defined  above,  these  columns  are  the  columns  of  T 
associated  with  the  zero  eigenvalues  of  S.)  Thus,  Equation  (22)  provides 
a  solution  to  Equation  (14)  for  the  n  y ,  in  terms  of  the  n  independent  (the 
last  p)  elements  of  z.  Moreover,  Equation  (22)  is  an  "algorithmic" 
expression  in  that  standard  numerical  procedures  exist  for  matrix  diagonal- 
ization,  eigenvalue  determination,  and  hence,  for  the  numerical  evaluation 
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of  che  nxp  matrix  T. 


In  index  notation.  Equation  (22)  may  be  written  as: 

y-  =  t.r2c  '  =  l ..... r, ;  r  =  1 . p  (23) 

where  the  t  may  be  thought  of  as  components  of  the  column  eigenvectors 

t  in  n  dimensional  space. 

Finally,  the  formal  statement  of  Equation  (22)  constitutes  the  "zero 
eigenvalues  theorem"  [23] . 


GOVERN IMG  EOUA 1  LONS 


The  procedures  outlined  above  can  be  used  to  systematically  formulate 
the  solution  to  the  multibody  system  dynamics  equations  (9)  subject  to  the 
constraint  equations  (13).  To  develop  this  formulation,  consider  again  the 
partial  velocity  and  partial  angular  velocity  vectors  discussed  above. 

From  Equations  (1)  and  (23)  the  velocity  of  in  R  and  the  angular  velocity 
of  in  R  may  be  expressed  in  the  form: 


v,  =  v,  t  ,  z  n  and  . ,  =  ....  t  .  on 
~k  k.n  ..  r  r~m  'k  k-m  r  r~m 


(24) 


where  the  z^_  (r=l,...,p)  may  be  viewed  as  new  generalized  coordinate 

derivatives.  The  partial  velocity  of  G^  in  R  and  the  partial  angular 

velocitv  of  3,  in  R,  with  respect  to  z  ,  then  become: 
k  r 


Hence , 


V’Zr 


v,  t  ,  n 
k>m  '.r~m 


and 


/  3  z  =  u  t ,  n 
~k  r  k- m  • r-m 


the  generalized  active  and  inertia  forces  of  Equations 


(25) 

(7)  and 


*  a 
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(8)  become: 


"  =  •’  t  F  c  M 

r  k;m  ir  1m  kirn  l r  km 


F*  =  v,  .  t ,  F*  +  o,  ,  t ,  M* 
r  k.a  ir  km  k’.m  ir  km 


Then,  from  Lagrange's  form  of  d'Alembert's  principle,  the  governing 
equations  (6)  become: 


F*  =  0  r  =  1 , . .  .  ,  p 

y  ' 


or,  more  specifically; 


:  •••.•=  f  .  t . 

car-  ■  r 


r  =  i .....  p 


where  a,  and  f l  are  given  by  Equations  (10)  and  (11). 

L  CJ  ^ 

Equations  (29)  together  with  the  constraint  equations  (13)  constitute 
the  system  of  equations  to  be  solved.  A  numerical  procedure  for  their 
solution  can  be  formulated  as  follows:  Consider  the  general  case  where  p=n-m. 
Then,  by  differentiating,  the  constraint  equations  (13)  become: 


b  V  “  —  b  V 

i  •  ■  ■  i  ' 


i  =  1 , .  .  .  ,  m 


Equations  (29)  and  (30)  form  a  total  of  n  equations  for  the  2n  unknowns 

y,  and  x,.  Hence,  there  needs  to  be  annexed  to  these  equations  the  expressions: 


for  the  consistent  numerical  formulation  of  the  governing  equations.  (If 
the  v  are  chosen  such  that  the  x  do  not  exist,  as  mentioned  earlier. 
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chen  Equations  (111  must  be  replaced  by  analogous  expressions  relating 
v  to  other  variables  (such  as  Euler  parameters  [10J  vhicn  define  the 
relative  orientations  of  the  bodies.) 

The  balance  of  the  numerical  formulation  of  the  solution  of  Equations 
(29),  (30),  and  (31)  is  now  routine:  It  is  perhaps  most  convenientlv 
expressed  in  matrix  notation.  To  this  end,  let  C  be  the  nxn  matrix 
containing  the  coefficients  of  y,  in  Equations  (29)  and  (30).  Then,  in 
partitioned  form  C  is: 


where  x  and  y  are  the  column  matrices  with  elements  x  and  v,  ( .’  -  l,...,u) 
respectively. 


14 


example 


For  a  simple  example  illustrating  some  ol  these  ideas,  consider  the 
planar  triple  pendulum  shown  in  Figure  3.  The  three  rods  are  identical 
having  length  •,  and  there  are  f ric cionless  pins  at  the  joints:  0^,  0^, 
and  0^.  The  system  has  3  degrees  of  freedom  which  may  be  described  by 
the  orientation  angles  0  ,  90,  and  9^  shown  in  the  figure.  Using  Lagrange's 
form  of  d'Alembert's  principle,  the  equations  of  motion  of  the  system 

V-.  -  IS,;,  *  hi’ ,  +  6C9+3)-'1  +  (10  +  9C,  ^  6C3  +  3C,+  j>  , 

+  (3  +  3C3  +  .-J  =  -(g/O  (15Sl  +  9S1+2  +  3S1+2+3) 

-  (9S2  +  3S2+3)0“  +  (9S2  -  3S3> (e3  +  02)2 

+  (3S2+3+  3S3)(‘l  ‘  *2  +  V"  (35'; 

(10  *  9C„  +  60,  +  +  (10  +  60-,)%  +  (2  +  3C-)0„ 

)  -I4-  3  !  3  l  3  j 

=  -(g/.)(9S1+,  +  3S1+2+3)  -  ( 9S ^  +  3S2+3^ei 

-  3S3(01  +  9,)2  +  3S3(9x  +  92  +  93)2  (36) 

and 

(2  -  3C3  +  3CO+3)01  +  (2  +  3C3)92  +  2?3 

=  3^g/  ,^1  +  7+3  -  -5S^_j_3j3  -  3S 3 ( ^ ^  +  (3/) 

where  C.  =  cos6 . ,  C...  =  cos (9.  +  3.),  etc. 

1  i  i+j  1  j 


15 


A  "closed  loop"  or  constraint  mav  be  formed  by  fixing  the  end  point 
P  of  the  pendulum.  Hence,  let  P  be  fixed  at  a  point  having  coordinates 
(a,b)  relative  to  the  X-Y  coordinate  system  shown  in  Figure  3.  This 
constrained  system  has  only  one  degree  of  freedom.  Two  scalar  constraint 
equations  relating  the  coordinates  S^,  and  maY  be  obtained  from  the 

position  vector  equation: 


0,0„  +  0-0.,  +  0,P  +  P  0.  =  0 

12  23  3o  o_  i 

That  is,  considering  the  horizontal  ar.c 

equation  leads  to  the  equations: 


(38) 


uciits  or  tnis 


SL  ~  Sl+2  +  H 1-2+3  a/' 


(39) 


'1  1+2  1+2+3 


=  b/e, 


which,  upon  differentiation,  become: 


)0,  +  +  c^+2+3^3  ® 


(C1  +  Cl+2  +  Cl+2+3)0l  +  (Cl+2  +  Cl+2+3 


(40) 


(S1  +  Sl+2  +  Sl+2+3)0l  +  (Sl+2  +  Sl+2+3)02  +  Sl+2+303  =  ° 


Equations  (39)  and  (40)  represent  Equations  (12)  and  (13)  in  the  foregoing 
analysis. 

To  simplify  the  analysis,  let  P  be  fixed  on  the  X-axis  at  a  point  Pq 
a  distance  ?.  from  0^  The  system  then  takes  the  form  of  a  rhombic  linkage 
as  shown  in  Figure  4.  In  this  case,  a= l ,  b=0,  and  the  constraint  equations 
(39)  are  seen  to  be  satisfied  by  the  relations:  90  =  n / 2—0 ^  and  3^  =  i/2+n^. 
The  coefficient  matrix  B  of  Equations  (13)  and  (40)  may  then  be  expressed  as; 
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1°  ' 

•ci 

1 

-C1  ! 

B  = 

! 

i 

i  (i  - 

■  V 

i 

~si 

j 

(41) 

The  matrix  S  of  Equation 

(15) 

then 

becomes : 

i 

(1 

-  V 

-S1 

S  = 

(i  -  sx) 

2(1 

-  sl> 

(l  -  sL) 

(42) 

-si 

(1 

-  V 

1 

It  is  readily  seer,  chat  S  has  one  zero  eigenvalue  and  chat  Che  associated 
eigenvector  array  T  is: 


1 


(43) 


The  governing  equations  to  be  solved  may  now  be  obtained  using  Equations 
(29)  and  (30).  From  Equations  (35),  (36),  (37),  and  (40),  these  become: 

(16  +  9C?  +  3C.  +  6C2+3)ei  +  (2  +  9C,  +  3C3  +  3C,+3)02 

+  (2  +  3c2+3)e3  =  -(g/i)(i5s1  +  3s1+2+3)  -  3s2+3eJ 

+  (9S2  -  3S3)(01  +  02)2  +  3S2+3(01  +  §2  +  03)2  (44) 


"C1S2  " 


C1S3  =  W2  +  si9r3 


(45) 


and 


(1  -  S1)e,  -  s103  =  -c15192 


W3 


(46) 
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Equations  (44),  (45),  and  (46)  have  been  solved  numerically.  The 


results  :natch  the  results  or  the  numerical  solution  of  the  pendulum 
equation: 

'3  ^  +■  (Sg/SiOsin^  =  0  (47) 

which  is  the  governing  equation  of  the  rhombic  linkage  of  Figure  4. 


REDUCTION  OF  THE  NUMBER  OF  GOVERNING  EQUATIONS 

A  principal  step  in  the  foregoing  formulation  as  well  as  in  the  above 
examp Le  is  the  differentiation  of  the  constraint  equations  (See  Equations 
(36),  (45),  and  (46).).  Since  these  equations  are  then  combined  with  the 
reduced  set  of  dynamical  equations  and  the  resulting  system  is  integrated, 
a  question  which  arises  is:  Is  it  necessary  to  first  differentiate  and 
later  integrate  these  equations?  That  is,  could  some  computational 
efficiency  be  obtained  by  avoiding  the  differentiation-integration  steps? 
The  answer  to  these  questions  is  that  it  is  indeed  possible  to  integrate 
fewer  equations  and  thus  obtain  some  computational  advantages.  A  procedure 
demonstrating  this  is  outlined  in  the  following  paragraphs: 

Let  the  mxn  matrix  B  in  Equation  (14)  be  partitioned  into  two 
submatrices  as : 

B  =  [B  :  BJ  (48) 

a  '  b 

where  B  is  an  m x  (n-m)  arrav  and  B,  is  an  m x m  square  arrav.  This 
a  b 

partitioning  of  B  induces  a  partitioning  of  the  v  array  in  Equation  (30). 
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That  is,  if  Equation  (30)  is  written  in  the  matrix  form  as:  Bv  =  -By, 


then  3v  maw  be  exoressed  as: 


*  *  l‘.  IV  ri-  * 


where  v  is  a  column  arrav  with  n-m  elements  and  y,  is  a  column  arrav  with 
a  b 

m  elements.  If  B  ^  exists.  Equation  (49)  may  be  solved  for  y  in  the  form: 
b  b 

yb  =  "Bb1(Va  +  By>  (50) 

Similarly,  let  the  reduced  set  of  dynamical  equations  (29)  be  written 

in  the  matrix  form  Av  =  f  where  A  is  the  (n-m)xm  arrav  with  elements  c  a 

rq  a 

r=l,,,.,n-m;  n  and  f  is  the  column  arrav  with  elements  f  t, 

■.  '  r 

r=l,...,n-m.  Then  following  the  pattern  of  Equations  (48)  and  (49>,  let 
the  dynamical  equations  be  written  in  the  partitioned  form: 


tAa  :  ‘V 


i—  — 
* 

ya 

•  i 


where  A^  is  an  (n-m)x(n-m)  square  array  and  A^  is  an  (n-m)  xm  array. 

By  substituting  from  Equations  (50)  and  by  matrix  block  multiplication 
Equation  (51)  becomes: 


A  y  +  A,  [-B,  (B  v  +  Bv)]  =  f 
a  a  b  baa  - 


By  rearranging  the  terms,  this  equation  may  be  written  as: 


Aya  =  f 


where  A  is  the  (n-m)  x  (n-m)  square  array  defined  as: 
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(54) 


and  f  is  the  (n-m)  element  column  array  defined  as: 

f  =  f  +  A^B^Bv  (55) 

b  b 

Equation  (53)  is  equivalent  to  a  system  of  (n-m)  scalar  differential 
equations  containing  2n  unknowns:  y^  and  x^ ,  £=l,...,n.  Hence,  there 
needs  to  be  annexed  to  this  system  an  additional  n-Hn  scalar  equations. 
Equation  (31)  provides  n  of  these  equations.  In  matrix  form  they  may  be 
written  as: 


x  = 


y 


(56) 


The  final  m  equations  may  be  obtained  from  Equations  (50).  However,  unlike 
Equations  (53)  and  (56) ,  Equations  (50)  are  algebraic  equations  and  do  not 
need  to  be  integrated.  That  is,  the  system  of  2n  equations  of  Equations 
(501,  (53),  and  (56)  contain  2n-m  differential  equations  and  m  algebraic 
equations  for  the  2n  unknowns  v„  and  xf ,  1*1,..., n.  This  is  a  reduction  of 
m  differential  equations  from  the  previous  system  of  Equations  (29),  (30), 
and  (31). 


DISCUSSION 

k. 

! 

•»  | 

At  this  point  there  are  several  comments  and  observations  which  might 

be  helpful.  First,  in  the  procedure  of  the  zero  eigenvalues  theorem,  the 

t 

m  constraint  equations  are  solved  for  the  n  y  ^  in  terms  of  n-m  new  variables 

z  .  Interestingly,  in  the  subsequent  formulation  of  governing  equations, 

%  r 

these  new  variables  z  do  not  appear.  Indeed,  it  is  onlv  the  coefficients 
■i 

«  ; 

- 
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t,  of  the  z  which  are  used.  As  noted  earlier,  these  coefficients  are  the 
components  in  n-aimensional  space  of  the  eigenvectors  t  associated  with  trie 
zero  eigenvalues  of  S.  However,  in  this  context,  since  the  corresponding 
eigenvalues  are  zero,  St^  is  zero  and  the  eigenvectors  tf  are  thus 
"orthogonal"  to  the  rows  of  S.  This  in  turn  means  that  these  eigenvectors 
are  orthogonal  to  the  rows  of  the  constraint  matrix  B.  (This  conclusion 
was  also  reached  in  an  earlier  analysis  of  constraint  equations  in 
n-dimensional  space  [27].)  Hence,  let  the  rows  of  B  be  thought  of  as 
"constraint  vectors"  in  n-di~e--.sion.ii  snace.  Then,  since  the  t  are  used 

■-V 

to  form  the  new  partial  velocitv  and  partial  angular  velocity  vectors,  the 
physical  system  can  be  considered  to  be  constrained  to  move,  in  n- dimension.-, 
snace,  in  directions  orthogonal  to  these  constraint  vectors  —  that  is. 
in  directions  defined  by  the  eigenvectors  t^. 

Next,  Lagrange's  form  of  d'Alembert’s  principle  is  an  ideally  suited 
method  for  formulating  the  dynamical  equations  when  there  are  accompanine 
constraint  equations.  Indeed,  the  governing  differential  equations  may  be 
developed  by  simply  contracting  the  dynamical  equations  obtained,  via  the 
principle,  by  using  the  t.,_  arrav  obtained  from  the  zero  eigenvalues  theorem 
This  procedure  is  seen  to  be  successful  since  the  generalized  forces  are 
linear,  homogeneous  functions  of  the  partial  velocity  and  angular  velocity 
vectors,  which  in  turn,  are  coefficients  of  the  generalized  coordinate 
derivatives  (in  the  velocitv  and  angular  velocity  vectors'!.  Therefore,  a 
modification  of  these  derivatives  directly  changes  these  vectors  and  hence, 
also  the  generalized  forces.  This  means  that  the  modification  procedure 
tor  the  generalized  coordinate  derivatives,  as  developed  by  the  zero  eigen¬ 
values  theorem,  may  be  directly  applied  to  the  dynamical  equations  them- 
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selves.  Also,  due  Co  these  arguments,  it  is  seen  that  ciiis  procedure  would 
not  necessarily  be  successful  if  the  dvnamical  euuations  were  obtained  by 
some  other  method.  (Additional  discussion  of  the  merits  of  Iterance's  form 
of  d'Alembert's  principle  may  be  found  in  references  [10,11,20,21,11].; 

Finally,  the  procedure  developed  herein  is  deemed  to  be  well  suited 
for  the  automated  development  of  the  governing  equations.  Numerical 
algorithms  are  currently  being  prepared  to  be  incorporated  into  the  computer 
codes  discussed  in  [10,11,12].  Additional  information  on  this  may  be 
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